#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <gmp.h>
#include <time.h>

#define SIZE_32 4
#define WORDS_32_BITS 32
#define NAILS 0
#define BASE 2
#define ORDER 1
#define ENDIAN 0
#define POWER_32 4294967296LLU

/*
 * Implementation of algorithm 1.4 of the text book.
 * Parameter :  a,b large intergers represented in 2^32 base array of words
 * 				"a" needs to be the larger or equal in size (no. of words) to "b"
 * 			    sa - size of a
 * 				sb - size of b
 * 				c - return value
 * Output : size of c
 * 
 * Author : Thilina Gunarathne
 * Indiana University CSCI B503 class project
 */
static size_t Add(unsigned int *a, unsigned int *b, unsigned int *c, size_t sa,
		size_t sb) {

	int i;
	unsigned long long carry=0;
	unsigned long long s, a_num, b_num;

	// do the addition up to the size of the smaller number
	for (i = 0; i < sb; ++i) {
		a_num = a[i];
		b_num = b[i];
		// All the variables needs to type "long long" for the result to be "long long"
		s = a_num+b_num+carry;
		/* Hmmm... Following proved to be faster than bit shifting */
		c[i] = s;
		carry = s/POWER_32;
	}
	
	// Do the addition of carry for the rest of larger number
	while (i < sa && carry) {
		a_num = a[i];
		s = a[i]+carry;
		c[i++] = s;
		carry = s/POWER_32;
	}
	// Just copy the rest of the "a" to the output
	while (i<sa) {
		c[i++]=a[i];
	}

	// If a carry remains then the size is one larger than "sa".
	if (carry) {
		c[sa] = carry;
		return sa+1;
	} else {
		return sa;
	}
}

/*
 * Implementation of Subtract algorithm derived from the addition algorithm
 * 
 * We are substracting "b" & "c" from "a".
 * Numbers needs to be such that sa>sb>sc  
 */
static size_t Substract(unsigned int *a, unsigned int *b, unsigned int *d,
		unsigned int *c, size_t sa, size_t sb, size_t sd) {
	int i;
	long long sum=0, a_num, b_num, d_num;
	u_int32_t sum1, sum2, borrow =0;

	// substract b&c from "a" untill the size of c is reached
	for (i = 0; i < sd; ++i) {
		a_num = a[i];
		b_num = b[i];
		d_num = d[i];
		sum = a_num -b_num-d_num+(sum>>32);
		c[i] = sum;
	}

	// substract rest of "b" from "a"
	for (i = sd; i < sb; ++i) {
		a_num = a[i];
		b_num = b[i];
		sum = a_num - b_num +(sum>>32);
		c[i] = sum;
	}

	// substract from rest of "a" if there remains a borrow
	while (i < sa && sum<0) {
		a_num= a[i];
		sum = a_num +(sum>>32);
		c[i++] = sum;
	}
	// Copy the remains of the "a" to the output variable
	while (i<sa) {
		c[i++]=a[i];
	}
	return sa;
}

static size_t Fast_Multiply(unsigned int *a, unsigned int *b, unsigned int *c,
		size_t sa, size_t sb) {
	unsigned int *u1, *u2, *v1, *v2, *w1, *w2, *w3, *w4, *d, *t1, *t2;

	size_t su1;
	size_t sv1;
	size_t st1;
	size_t st2;
	size_t sw1;
	size_t sw2;
	size_t sw3;
	size_t sw4;
	size_t div;

	if (sa > 20&& sb>20) {
		if (sa<sb) {
			div = sa/2;
		} else {
			div = sb/2;
		}
		
		// Following break up the input number using pointers to the middle of the array
		u2 = a;
		u1 = a+ (div);
		su1 = (sa-div);

		v2 = b;
		v1 = b+ (div);
		sv1 = (sb-div);

		t1 = (unsigned int *)calloc(sizeof(int), su1+sv1+2);
		t2 = t1+su1+1;

		/* Step 1 of Algo 5.2 */
		// U1 is either large or equal in size to U2 according to the above division logic
		// Hence U1 as the first argument
		st1 = Add(u1, u2, t1, su1, div);

		/* V1 is the large number according to the above division logic */
		/* Step 2 of Algo 5.2 */
		st2 = Add(v1, v2, t2, sv1, div);

		/* 
		 * Step 3 of Algo 5.2   W3 --> T1 * T2 
		 */
		w3 = (unsigned int *)calloc(sizeof(int), st1+st2);
		sw3 = Fast_Multiply(t1, t2, w3, st1, st2);

		/* 
		 * Step 4 of Algo 5.2   W2 --> U1 * V1 
		 */
		w2 = (unsigned int *)calloc(sizeof(int), su1+sv1);
		sw2= Fast_Multiply(u1, v1, w2, su1, sv1);

		// Use the output array to construct the output
		// This avoids unnecesary copying at the end of program 
		w4=c;
		
		/* 
		 * Step 5 of Algo 5.2   W4 --> U2 * V2 
		 */
		sw4 = Fast_Multiply(u2, v2, w4, div, div);

		/* 
		 * Step 6 of Algo 5.2   W3 --> W3 - W2 - W4 
		 */
		//Making sure we pass the smaller size number as the 3rd parameter
		if (sw2>sw4) {
			sw3 = Substract(w3, w2, w4, w3, sw3, sw2, sw4);
		} else {
			sw3 = Substract(w3, w4, w2, w3, sw3, sw4, sw2);
		}

		/* 
		 * Step 7 of Algo 5.2   C <-- [W4/2^n] & W4 <-- W4 mod 2^n 
		 */
		d = w4+ (div);

		/* 
		 * Step 8 of Algo 5.2   W3 <-- W3+C  , C <-- W3 / 2^n & W3 <-- W3 mod 2^n 
		 */
		sw3 = Add(w3, d, c+(div), sw3, (sw4 - div));

		/* 
		 * Step 9 of Algo 5.2   W2 <-- W2+C  , W1 <-- W2 / 2^n & W2 <-- W2 mod 2^n 
		 */
		sw2 = Add(w2, c+(2*div), c+(2*div), sw2, (sw3 - div));

		free(w2);
		free(w3);
		free(t1);
		//calculate the size of the return number
		return (3* (div)) + (sw2 - div);
		
	} else {
		
		/*
		 * Multiply32 :Grade school algorithm 
		 * Algorithm 1.8 in the text book
		 */
		unsigned long long carry, y_num, z_num, x_num;
		unsigned long long p=0;

		int i, j;
		carry = 0;
		//Multiply the "b" by the first digit of "a". 
		//No need to do the addition of old value as the result is empty 
		if (a[0] !=0) {
			x_num = a[0];
			for (j = 0; j < sb; ++j) {
				p = x_num * b[j]+ carry;
				c[j] = p;
				carry = p/POWER_32;
			}
			c[sb] = carry;
		}

		for (i = 1; i < sa; ++i) {
			carry = 0;
			// Just pass on if a[i] is 0.
			if (a[i] !=0) {
				x_num = a[i];
				for (j = 0; j < sb; ++j) {
					z_num = c[i+j];
					p = x_num * b[j]+ z_num + carry;
					c[i+j] = p;
					carry = p/POWER_32;
				}
				c[i+sb] = carry;
			}
		}

		if (carry=0) {
			return i+j;
		} else {
			return i+sb;
		}
	}
}

// Stub function we got from the given framework
void Product32(u_int32_t *a, u_int32_t *b, u_int32_t *c, size_t sa, size_t sb,
		size_t *sc) {
	
	/******************* ADDITION TEST ************/
	//    	z_size[j] = Add(x[j], y[j], z[j], x_size[j], y_size[j]);
	/******************* End Addition Test *********/

	/******************* SUBSTRACT TEST ************/
	//    	z_size[j] = Substract(z[j],x[j], y[j], z[j],z_size[j], x_size[j], y_size[j]);
	/******************* End Addition Test *********/
	
	*sc = Fast_Multiply(a, b, c, sa, sb);
}
